Photoelectron spectra of anionic sodium clusters from time-dependent 
density-functional theory in real-time 
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We calculate the excitation energies of small neutral sodium clusters in the framework of time- 
dependent density-functional theory. In the presented calculations, we extract these energies from 
the power spectra of the dipole and quadrupole signals that result from a real-time and real-space 
propagation. For comparison with measured photoelectron spectra, we use the ionic configurations 
of the corresponding single-charged anions. Our calculations clearly improve upon earlier results for 
photoelectron spectra obtained from static Kohn-Sham eigenvalues. 

PACS numbers: 36.40.Cg,31.15.Ew,33.60.-q,61.46.-w 



I. INTRODUCTION 

Since more than one hundred years, photoelectron 
spectroscopy plays an important role in physics. Ein- 
stein's explanation of the photo-effect, probably the most 
well-known experiment in the field, was a crucial step in 
the development of quantum mechanics. Whereas the 
photo-effect revolutionized the understanding of light, 
the main aim of modern photoelectron spectroscopy is 
to understand electronic and ionic structures from solids 
down to single atoms. Especially in the context of 
nanoscale materials photoelectron spectroscopy is one of 
the most important experimental tools since it is almost 
the only method that provides access to the electronic 
and ionic structures of these materials. The direct obser- 
vation of the electronic shell structure in sodium clusters 1 
is just one example for the power of the method. Another 
application is the determination of the ionic structure of, 
e.g., clusters. Since the electronic structure, and thus the 
photoelectron spectrum (PES), depends on the ionic con- 
figuration, comparing the measured PES with the results 
from first-principle calculations allows the identification 
of the ionic structure. This interplay between theory and 
experiment has already been used successfully in many 
cases^^^OO. 

Clearly, the just mentioned method can only work if 
reliable calculations for the system of interest can be per- 
formed. Since most of the measured systems consist of 
many electrons, density-functional theory (DFT) is an 
especially well-suited tool due to its low numerical costs. 
Unfortunately, evaluating the PES from a Kohn-Sham 
(KS) DFT calculation is not an easy task since only the 
highest occupied KS eigenvalue has a rigorous connection 
to the PES: it is equal to the ionization potentia l 11 ! 12 ' 13 . 
Thus, it yields the position of the first peak in the PES. 

The most common approach to obtain the other peaks 
in the PES from a DFT calculation is based on the den- 
sity of states of the occupied KS orbitals. In this ap- 
proach the KS eigenvalue spectrum of the 'mother' sys- 
tem, i.e., the system which still contains the photoelec- 
tron, is directly compared to the experimental PES. In 
many situations the resulting spectra reproduce the ex- 



perimental ones quite we ll2iMi5iSiLSi2ii2iil^. 

Another way to extract the information related to the 
PES from a DFT calculation is via the excitation ener- 
gies of the 'daughter' system, i.e., the system with one 
electron less. Since time-dependent DFT (TDDFT)i&ll, 
in principle, allows to calculate the excitation energies of 
a system exactly, TDDFT can be used to calculate the 
positions of the PES peaks accurately. This approach is 
followed in Ref. 0J and Ref.JlS The basic idea of this 
approach is also used in Ref. [20j, but in combination with 
configuration- interaction and not TDDFT calculations. 

We finally want to mention a third method how the 
PES can be obtained from a DFT calculation. In this 
approach, the time-dependent ionization process is simu- 
lated in real-time and the kinetic energy spectrum of the 
outgoing components of the KS Slater determinant are 
connected to the PES^.. 

In Ref. 0, the PES resulting from Na^, Naf , and 
Na^~ (among others) irradiated by an XeCl excimer laser 
(Hlu — 4.02 eV) were measured and compared with 
the corresponding KS eigenvalue spectra. Although the 
agreement between the theoretical results and the mea- 
sured PES was generally reasonable, a systematic dis- 
crepancy was found. Namely, the width of the theoret- 
ical spectrum, i.e., the difference between the energy of 
the energetically highest and lowest occupied KS eigen- 
value, was too large by about 0.2-0.4 eV. In Ref. [22| the 
reason for this discrepancy was examined. It was shown 
that technical aspects, e.g., the treatment of the pseu- 
dopotential, could not explain the differences. Further- 
more, it was demonstrated that using the exchange-only 
optimizcd-cffcctivc potential (OEP) 23 reduced the width 
of the KS spectrum but not to an extent that would bring 
the spectrum in agreement with experiment. In addition, 
it was also shown that using Slater's transition state con- 
cept could also not improve the theoretical results. Thus, 
the question arises whether a different method to extract 
the PES from a DFT calculation leads to a better agree- 
ment with the experiment. It is the aim of the present 
manuscript to answer this question by extracting the PES 
from the excitation energies of the corresponding 'daugh- 
ter' systems. 
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FIG. 1: Schematic view of two different approaches to cal- 
culate the PES. Left: The process is described as a strong 
excitation of the (N + l)-electron system. Right: The pho- 
toelectron has already been detected and the remaining N- 
electron system is left in an excited state. The link between 
the kinetic energy of the photoelectron and the energy of the 
excited state of the iV-electron system is provided by energy 
conservation. 

II. THEORETICAL BACKGROUND 

Before discussing the results let us sketch the theoreti- 
cal background of the presented calculations in more de- 
tail. Fig. Q] schematically shows two approaches how the 
peak positions in the PES can be calculated. On the left 
hand side, the process is described as an excitation pro- 
cess from the ground state to an energetically high-lying 
state with continuum contributions. Since KS eigenvalue 
differences are zeroth-order approximations to excitation 
energies 2 ^^, the KS DOS of the (./V + l)-electron system 
can be used to obtain an approximate PES. In addition 
to this argument, Chong et alM> have given well founded 
arguments that KS eigenvalues can be interpreted as ap- 
proximations to relaxed vertical ionization potentials. 

On the right hand side of Fig. [TJ the situation after the 
photoelectron has been detected is considered. In this 
case, the remaining system is left in an energetically low- 
lying excited state of the iV-electron system. To connect 
the excitation energies of this system to the PES, energy 
conservation is used. Before the photon is absorbed the 
total energy is given by Eq N+1 ^ + Hlu, where Eq N+1 ^ is the 
ground-state energy of the 'mother' system containing 
N + 1 electrons and fiu> is the photon energy. After the 
detection of the photoelectron the total energy is given 
by the kinetic energy of the photoelectron -Ekin and the 
energy of the remaining 'daughter' system. Since the 
total energy is conserved, it follows that 

-Ebindj = -Ekin — hu> = Eg N+1 ^ — Eq N — AEj N (1) 



must hold. Here, Eq ' is the ground-state energy of 
the 'daughter' system and AEj are its excitation 
energies^ 2 . For the first peak in the PES the kinetic 
energy of the photoelectron is maximal. In this case, the 
'daughter' system is in its ground-state, i.e., AE^ is 

zero and the peak position is at Eq N+1 ^ — Eq . 

To obtain the excitation energies from time-dependent 
DFT the full linear density-response function of the in- 
teracting system can be used. This function provides ac- 
cess to the excitation energies of the system since it has 
poles at these energies. The crucial observation now is 
that the interacting linear density-response function can 
be expressed in terms of the KS response function and 
the exchange-correlation (xc) kerne l 17 ' 26 ' 27 . Nowadays, 
most applications use the matrix equation of Casida^ 7 - to 
obtain these excitation energies. 

Alternatively, the excitation energies can be extracted 
from a spectral analysis of the time-dependent density 
coming from a real-time propagatio n 28 ' 29 ' 30 . In this ap- 
proach, the xc kernel is not needed, but instead, the time- 
dependent KS equations are solved without explicit lin- 
earization. To illustrate this approach, imagine we have 
created a time-dependent density n(r, t) of an interacting 
system by, e.g., a laser excitation. Assuming that the sys- 
tem is confined by the same time-independent potential 
before and after the laser pulse we can write the excited 
density in terms of the eigenstates 1^-) of the interacting 
system in the time-independent potential. It reads 

n{r,t) = (V>(t)|r#(t)> 

= J2 c *j c ^j\^k) exp(-t(E k - E^t/h) (2) 

j,k 

where Ej is the eigenvalue corresponding to \tfij) and n is 
the density operator. Assuming that the time-dependent 
state \ip(t)) is dominated by the ground state, i.e., cq » 
Cj we can write 

n(r,t) w \c \ 2 n (r) + 

J2j cScj(^o|^|^j)exp(-i(Ej - E )t/h) + c.c.(3) 

Here, n (r) is the ground-state density of the system. If 
we now calculate the Fourier transform of n(r, t) we will 
get peaks at the exact excitation energies of the system. 
Since time-dependent DFT in principle provides us with 
the exact time-dependent density, this is an easy method 
to obtain the excitation energies of the interacting system 
from a time-dependent DFT calculation. 

In a practical calculation, two problems must be solved 
to get the excitation energies from this scheme. First, one 
has to create a time-dependent density which is domi- 
nated by the ground-state density and, in addition, con- 
tains the excited states of interest. The second problem 
is how to extract the excitation energies from the time- 
dependent density in practice. Since the density in every 
space point at all times cannot be stored, a full Fourier 
transform of Eq. ([3]) giving n(r, ui) is not possible. To 
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overcome this problem several possibilities exist. One 
is to evaluate n(r,u>) only for some points in spaced, 
e.g., in the center of the cluster. A different method is 
to Fourier transform certain moments of the density dis- 
tribution. Typically, the dipole moment is used for this 
purpos o 28 ! 29 . Obviously, some excitations are filtered out 
by this procedure because the Fourier spectrum of the 
dipole moment only shows excitation energies of states 
which can be coupled to the ground state via the dipole 
operator. Non-dipolc active excitations can be taken into 
account by recording higher moments. In the following 
we will use the timc-dcpcndcnt dipole and quadrupole 
moments obtained from a real-time propagation to ob- 
tain excitation energies for the systems of interest. 



III. TECHNICAL ASPECTS 

For the ionic ground-state configurations of the 
'mother' systems we used optimized structures obtained 
with the PARSEC 33 program package. The generalized- 
gradient approximation of Perdew et al. (PBE) 3 — was 
employed for the geometry optimization. The ionic cores 
were treated consistently with norm conserving non-local 
pseudopotentials 35 . 

The time-dependent KS equations were solved on a 
real-space grid in real time. For this we implemented 
the necessary algorithms into a modified version of the 
PARSEC code. In detail, we implemented a fourth-order 
Taylor approximation to the propagator in combination 
with a higher-order finite-difference formula for the ki- 
netic part of the KS Hamiltonian. For the calculations, 
we used a time step of 0.003 fs and the total propa- 
gation time was 75 fs. The ionic cores were again de- 
scribed by norm conserving non-local pseudopotentials. 
Furthermore, the ionic structures were fixed during the 
propagation. The grid spacing was 0.7 ao and the grid 
radius varied between 20 and 23 ao depending on the 
system. The time-dependent density was created by ap- 
plying a boost exp(zr • Pboost/^) to the ground-state 
KS orbitals. The total excitation energy of the sys- 
tem was -Ec xcit = 1-0 x 1 0~ 5 eV, i.e., a boost strength 
bboost | = \/2m e E CXC it / N was applied to each KS orbital 
(with m e being the electron mass and N the number of 
electrons). In addition, the calculations were repeated 
with a boost strength reduced by a factor of 1.0 x 10~ 2 . 
Using these two small boost strengths allowed us to check 
whether the created time-dependent density was domi- 
nated by the ground-state density (see below). 

Instead of applying the same boost vector Pboost to all 
KS orbitals, and thus creating a coherent velocity field, 
we varied the boost direction for different KS orbitals. 
This is necessary since applying the same boost direction 
to all KS orbitals corresponds to first order in Pboost to 
a dipole excitation of the system, i.e., from the result- 
ing time-dependent density it is only possible to retrieve 
the excitation energies of 'dipole-active' states. By ap- 
plying different boost directions to different KS orbitals 



we modeled a general excitation mechanism creating a 
time-dependent density containing excited states with 
different symmetry properties. In detail, we randomly 
chose a boost direction (no symmetry axis of the con- 
sidered cluster) for the first orbital and then chose our 
coordinate system such that this direction was the first 
diagonal (for the remaining rotational degree of freedom 
a random angle was chosen). After this we boosted the 
second orbital in the opposite direction of the first boost. 
The third orbital was then boosted in the direction of 
the second diagonal of the chosen coordinate system, the 
forth again in the opposite direction and so on. For Nag, 
the ninth orbital was boosted again in the same direction 
as the first orbital. Since the only purpose of this pro- 
cedure was to create a time-dependent density without 
any particular symmetries, we do not consider the rela- 
tive orientation of the cluster with respect to the boost 
directions to be of special importance. 

Finally, we used the time-dependent local-density ap- 
proximation (TDLDA) 36 for the xc potential for the 
propagation. Since the linear response of the homoge- 
neous electron gas is the same in this approximation and 
in the PBE functional, the differences in the resulting 
excitation energies can be expected to be small in the 
low-energy regime. Unfortunately, the exact-exchange 
orbital functional cannot be used for comparative calcu- 
lations since it requires a solution of the time-dependent 
optimized-effective potential equation^! and no method 
for this exists in real-time at the moment 38 . 



IV. RESULTS AND DISCUSSION 

A. Results for Na^ 

Fig. [2] shows the dipole power spectra of Na3 resulting 
from two boost strengths differing by a factor of 10 2 . The 
dipole power spectrum is given by 

:= 5>i( w )| a (4) 

3=1 

with dj(oj) being the Fourier transform of the dipole mo- 
ment 

dj(t) = J Xjn{r,t) d 3 r (5) 

where x\ corresponds to the Cartesian coordinate x, xi 
to y, and X3 to z. For small momentum boosts first- 
order perturbation theory predicts a linear dependence 
of the expansion coefficients Cj in Eq. {2| on the boost 
strength. As a consequence, reducing the boost strength 
by a factor of c suppresses peaks corresponding to energy 
differences between two excited eigenstates by a factor of 
c 4 in the power spectrum. Since peaks corresponding 
to transitions between the ground state and an excited 
eigenstate are only suppressed by a factor of c 2 , chang- 
ing the boost strength allows one to distinguish between 
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FIG. 2: (Color online) Dipole power spectrum of Na3 resulting 
from an incoherent boost excitation. The result obtained from 
a total excitation of 1 x 10" 
whereas the label 'weak excit.' 

by a factor of 1 x 10 -2 . Clearly, the dipole power spectrum 
scales quadratically with the boost strength indicating that 
the peak positions correspond to excitation energies between 
the ground state and excited states. 



5 eV is labeled 'strong excit.' 
corresponds to a boost reduced 
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FIG. 3: (Color online) Sum of the absolute square of the 
Fourier-transformed components of the quadrupole tensor re- 
sulting from the same excitations as in Fig. [2] In contrast to 
the dipole power spectrum some peaks vanish with reduced 
boost strength indicating that they correspond to energy dif- 
ferences between excited states. 



these two kinds of excitations. As one can see in Fig. [2j 
the results for the two boost strengths are almost iden- 
tical except for the predicted factor of 10 4 . Thus, we 
conclude that all the peak positions in the dipole power 
spectrum of Fig. [2] correspond to energy differences be- 
tween the ground state energy and the energy eigenvalues 
of the excited eigenstates. 

The situation is different for the power spectrum re- 
sulting from the quadrupole moments. In Fig. [3] we plot 



i=l 
3>i 



(6) 



for the same two excitation boosts. In this equation, 
qij (u) is the Fourier transform of the quadrupole moment 



qij(t) 



n(r,t)(3xiXj - r 2 Sij) d 3 r , (7) 



and the sum only runs over the independent components 
of the quadrupole tensor. Clearly, the quadrupole spectra 
for the different excitation strengths differ considerably. 
For instance, the three large peaks at around 0.3, 1.1, 
and 1.7 eV vanish almost completely. Thus, we conclude 
that they belong to transition energies between different 
excited states. Indeed, one can see that these energies 
are exactly equal to the energy differences between the 
first excited state and the other excited states from the 
dipole spectrum. 

The reason why there are no peaks at these energies 
in the dipole spectrum can easily be understood if one 
takes the geometry of Na3 into account. Since Na3 has a 
linear ionic configuration, the ground state has even par- 
ity. Thus, the dipole spectrum only shows excited states 
with odd parity. Since two states with odd parity cannot 
be coupled by the dipole operator, transitions between 
these states do not show up in the dipole spectrum. 

After the identification of the true excitation energies 
we can now compare the results with the measured PES. 
In Fig. Q]the excitation energies of Na3, the KS DOS of 
Na^ and the measured PES (of Na^) are plotted. The 
positions of the occupied KS eigenvalues are indicated 
by arrows, long bars indicate excitation energies from 
the dipole spectrum and shorter bars excitation energies 
from the quadrupole moments. In addition, excitation 
energies leading to peaks below the strongest bound ex- 
perimental peak are reduced in their overall height. For 
better comparison, the KS DOS and the excitation spec- 
trum are both rigidly shifted in such a way that the most 
weakly bound peak coincides with the experimental one. 

As the upper part of Fig. 0] shows the peak positions 
that one obtains from the KS DOS are close to the ex- 
perimental peak positions. Unlike in the case of larger 
Na clusters, the width is slightly smaller than the en- 
ergy difference between the two large experimental peaks 
but it is still reasonable. However, since there are only 
two occupied KS orbitals in Na^ the KS DOS picture 
fails completely to describe the higher lying peaks in the 
measured spectrum. 

As one expects from Eq. (Q} the PES obtained from the 
excitation energies shows a much richer structure than 
the KS DOS. One striking feature for instance is the sec- 
ond excitation around —2.0 eV. It seems that the energy 
difference between this peak and the one at —1.7 eV is 
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FIG. 4: (Color online) Measured PES of Na^ ('Exp.') and 
theoretical PES obtained from the excitation energies of Na3 . 
Excitation energies from the dipole spectrum are labeled 'Dip. 
excit.' whereas 'Quad, excit.' labels excitation energies de- 
duced from the quadrupole moments. Arrows indicate the re- 
sult obtained from the KS DOS. Upper part: results obtained 
from the ionic ground-state configuration at zero tempera- 
ture. Lower part: results obtained from an ionic configuration 
with a larger bond length to simulate a higher temperature. 
For most peaks the agreement with the experimental PES is 
clearly improved. 



too large in the calculation and that they are merged to 
one peak in the experimental PES. However, in general, 
the dynamically calculated excitation energies and the 
energies obtained from the experimental PES are close 
to each other even for the stronger bound peaks. To see 
if the remaining discrepancy can be further reduced by 
taking temperature effects into account we have repeated 
our calculations with a larger bond length. Due to the 
net negative charge of the cluster one can expect that 
other geometry changes, e.g., bending, only play a minor 
role in the case of Na^~ . We have used a new bond length 
of approx. 6.8 instead of 6.5 do- This new value for the 
bond length I of the cluster has been obtained from an 
estimate for the thermal expansion at T = 300 K. It is 
based on the formula (3 = j-^p for the linear thermal 
expansion coefficient (3 which we have roughly estimated 
by P » 2 /3 b uik (see Ref. 39), where (3 bulk is the bulk value 
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FIG. 5: (Color online) Same as in Fig. Q] but for Na^ . Al- 
though both, the KS DOS and the PES from the excita- 
tion energies, describe the measured PES acceptable the large 
peak at — 2.2 eV is much better reproduced by the excitation 
energies from the 'daughter' system. The experimental data 
is taken from Ref. [y. 



for crystalline sodium at room temperature. 

The result can be seen in the lower part of Fig. 2] 
For most peaks one can observe a small shift towards 
lower absolute binding energies. Except for the peak at 
— 1.7 eV the agreement between the experimental and 
theoretical spectrum is slightly improved by the increased 
bond length. Especially the broader peak at around —2.6 
eV is nicely reproduced in this case. All in all, both 
calculations show that for Na^~ the main advantage of 
the 'excitation picture' is the reproduction of the deeper 
bound structures in the PES. 



B. Results for Na 5 

Fig. shows the experimental PES of Na^~, the KS 
DOS, and the PES obtained from the excitation energies 
of Na5. The labeling is the same as in the corresponding 
previous figures. As for Na^~, the KS DOS is in accept- 
able agreement with the first large peaks, although the 
strongest bound large peak has a too negative binding en- 
ergy in the KS DOS. As one can see these peaks are also 
well described by the excitation energies of the 'daugh- 
ter' system with the additional advantage that the last 
peak at —2.2 eV is better reproduced. In this approach 
it consists of four close-lying excitations. 

Beyond the peak at —2.2 eV the comparison with the 
experimental measurement is difficult since no clear peak 
structures can be observed. Perhaps the accumulation of 
excited states around —3.3 eV can be associated with the 
measured peak in this region, but for the reasons given 
below, one has to be very cautious in making comparisons 
in this part of the spectrum. 

As one can see from the results for Na^T and Nag dis- 
cussed below the problem of comparing the deeper lying 
part of the measured PES with calculated excitation en- 
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ergies is not specific to Na^~. In general, the density 
of excited states grows with the excitation energy, i.e., 
more and more states appear in the theoretical calcula- 
tion. On the other hand, from the point of view of first- 
order perturbation theory the PES depends not only on 
the positions of the excited states but also on the matrix 
element of the perturbing operator D between the ini- 
tial and the final state. Taking the ground state of the 
'mother' system for the initial state and a product state 
consisting of one photoclcctron with momentum k for 
the final state one obtains matrix elements of the form 
(k, ipj \ D \iJ)q N+1 '}. It is intuitively clear that these 
matrix elements are much larger for low-lying states than 
for energetically high-lying ones which in an independent- 
particle picture would correspond to removing one parti- 
cle and exciting a second one above the Fermi level. Espe- 
cially, in the case of truly independent particles this pro- 
cess cannot happen if the perturbing operator is a one- 
particle operator like the dipole operator. Thus, many 
energetically high-lying eigenstates of the 'daughter' sys- 
tem are hardly or even not at all excited in the experi- 
ment. Since the mentioned matrix elements depend on 
the interacting many-particle wavefunctions, calculating 
these exactly is close to being impossible. Especially, 
retrieving these matrix elements from a TDDFT propa- 
gation of the 'daughter' system is not trivial because the 
propagation only provides information about matrix ele- 
ments between excited states and the TV-particle ground 
state and not the (N + l)-particle ground state. 

However, as the presented calculations show, the ma- 
trix elements do not play a very important role in the 
part of the spectrum that we are mainly interested in. 
Nevertheless, the calculations also clearly indicate that 
one has to consider them if the deeper lying parts of the 
spectrum are of interest. A possible method how this can 
be done in a TDDFT calculation can be found at the end 
of Sec. El 



C. Results for Na 7 

The results for Naf are shown in Fig. [5J As said pre- 
viously, in the region below —2.5 eV, it is difficult to 
compare theory and experiment due to the great num- 
ber of close lying transitions. As in Na^, the KS DOS 
describes the strongest bound large peak worst. In this 
case it is already off by 0.4 eV. In contrast, the peak po- 
sition obtained from the TDLDA excitation energies is 
considerably closer to the experimental peak. It is only 
off by 0.1 eV. Thus, the overestimation of the spectrum's 
width by the KS DOS&22 is not observed in the result 
obtained from the TDLDA calculation. The remaining 
difference of 0.1 eV between the width of the theoreti- 
cal and the experimental spectrum can be easily caused 
by technical aspects like the employed pseudopotential 
and xc potential^. In addition, thermal effects like bond 
elongation and structural isomerization can shift the ob- 
tained width by 0.1 eV— i 39 i 40 . Considering that the ex- 
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FIG. 6: (Color online) Same as in Fig. [4] but for Naf . Espe- 
cially, the peak at —2.4 eV is much more accurately described 
by the excitation energies than by the KS DOS. In the weakly 
bound region thermal effects play an significant role in the 
case of Naf . This explains the rather poor agreement be- 
tween the theoretical values calculated at zero temperature 
and the measured curve between —1.3 and —1.7 eV. The ex- 
perimental data is taken from Ref. 0. 



perimental PES were obtained from clusters with a tem- 
perature of around 250-300 K, the difference between the 
theoretical result at zero temperature and the experimen- 
tal result is hardly surprising. At these temperatures the 
larger anionic sodium clusters behave liquidlike^. Conse- 
quently, many different ionic configurations are present 
in the experiment and show up in the measured PES. 

This aspect must also been kept in mind if the theo- 
retical and experimental results are compared in the re- 
gion between —1.3 and —1.7 eV. In this region both the 
zero temperature KS DOS result and the zero tempera- 
ture result from the excitation energies do not describe 
the measured PES very accurately. Especially, the ex- 
citation peak at —1.45 eV does not fit very well. How- 
ever, from Ref. it is known that the agreement be- 
tween the experimental and the KS DOS result in this 
energy region is significantly improved if different ionic 
structures are taken into account via Born-Oppcnheimer 
Langevin molecular-dynamics^. Therefore, one can ex- 
pect that the agreement between the experimental and 
the TDLDA result is also improved if different ionic struc- 
tures are taken into account. Due to the more compli- 
cated structures and the growing number of isomers the 
inclusion of the temperature influence on the ionic struc- 
tures of larger clusters is much more involved than in the 
case of Na^ . Additionally, combining Born-Oppcnheimer 
Langevin molecular-dynamics with the calculation of ex- 
citation energies is substantially more expensive than 
combining such a molecular dynamics scheme with a KS 
DOS calculation. Thus, including thermal effects in the 
present study is beyond the scope of the present work 
and is a future project. 
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FIG. 7: (Color online) Same as in Fig. [4] but for Nag . As 
in Naf , especially the stronger bound part of the spectrum 
is described more accurately by the calculated excitation en- 
ergies than by the KS DOS. The experimental data is taken 
from Ref. d 



D. Results for Na„ 



Finally, the theoretical results for Na^~ are compared in 
Fig. [7] with the measured PES. This cluster is the first one 
which has a clear peak in the range between the highest 
and lowest occupied KS eigenvalue which is completely 
absent in the KS DOS, i.e., the experimental PES shows 
six clear peaks whereas the KS DOS consists of only five 
peaks: the peak around —2.4 eV is completely missing 
in the KS DOS. In addition, the strongest bound peak 
in the KS DOS is off by 0.5 eV. In other words, the KS 
DOS result is inaccurate to the extent of being useless 
below —2.2 eV. As for Na^T , the splitting of the large peak 
around —1.8 eV is reproduced if different ionic structures 
are used£. 

In contrast to the KS DOS result, the PES obtained 
from the excitation energies is close to the measured 
curve over the whole range. Below the lowest lying peak 
at —2.7 eV the comparison is again difficult without 
knowing the matrix elements mentioned above. For the 
two peaks at —2.7 and —2.4 eV the theoretical values are 
off by 0.1 eV. Especially since Na^~, in contrast to Nay , is 
not a closed-shell cluster, one can expect that such energy 
differences can be easily caused by ionic structure modifi- 
cations induced by finite temperatures. As expected from 
Ref. the splitting of the peak at —1.8 eV is also not 
reproduced by the zero temperature TDLDA calculation. 
All in all, the experimental result in the weaker bound 
part of the spectrum is described equally well by the KS 
DOS and the excitation energies of the 'daughter' system. 
However, in the stronger bound part the time-dependent 
calculation yields a much more realistic description of the 
PES than the KS DOS. Since this emerges as a general 
observation for all systems studied in this manuscript, we 
discuss it on general grounds in the following section. 



V. SUMMARY AND CONCLUSION 

Using TDDFT we have calculated the excitation en- 
ergies of small neutral sodium clusters. The energies 
of the excited states were retrieved from the dipole and 
quadrupole moments of the time-dependent density via 
spectral analysis. The time-dependent density was cre- 
ated by an incoherent boost of all KS ground-state or- 
bitals and then propagated in real-space and real-time. 
To discriminate between true excited state energies and 
energies corresponding to energy differences between ex- 
cited states we did two calculations for all systems using 
different excitation energies. For comparison with mea- 
sured PES of the anionic clusters the excitation energies 
were calculated in the ionic configuration of the anions. 

In general, the PES for all clusters studied in this 
manuscript can be divided into three parts. The first 
part consists of large, 'weakly' bound peaks, the second 
of large, 'strongly' bound peaks, and finally a 'less struc- 
tured' region below the lowest large experimental peak. 
Except for Na^~ no comparisons between our theoretical 
and the experimental results can be made in the third 
region. As discussed in Subsection IIV Bl the main reason 
for this is the missing access to the transition matrix ele- 
ments between the ground and the excited states. Since 
the number of excited states can grow very rapidly, one 
can expect that the omission of the transition matrix ele- 
ments can cause severe problems if more complex systems 
are examined. A possible way to overcome this problem 
is by including the information of the matrix elements in 
the initial density, i.e., by creating an initial density that 
only includes the states which are really excited in the 
ionization process. Work in this direction is under way. 

In the middle part of the spectrum the results obtained 
from the TDLDA excitation energies are clearly superior 
to the results from the KS DOS. Especially, the position 
of the strongest bound large peak is much better repro- 
duced by the TDLDA calculation than by the KS DOS. 
Thus, using the TDLDA cures the main problem that 
plagues theoretical results obtained from the KS DOS 
for sodium clusters, namely the prediction of a signifi- 
cantly too large width of the spectrum. In addition, the 
PES from the TDLDA excitations can describe an ex- 
perimental peak in the PES of Na^~ which is completely 
missing in the KS DOS. The remaining differences be- 
tween the experimental and our theoretical results are all 
small enough to be explainable by technical details or the 
finite temperature (250-300 K) of the ionic structures in 
the experiment. In particular the finite temperature can 
be expected to be responsible for the difference since the 
considered clusters behave liquidlike at this temperature 
and thus, the measured PES result from many different 
ionic structures / isomers which differ from the theoreti- 
cal zero temperature ground-state structures used for the 
calculations. 

Finally, in the most weakly bound part of the spectrum 
we find that the TDLDA result and the one from the KS 
DOS are very similar. Since the KS DOS at finite tern- 
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perature is in very good agreement with the experimental 
result^, it is extremely likely that also the TDLDA exci- 
tation energies calculated from higher temperature ionic 
structures will describe the experimental PES very well 
in this region. 

Generally, our findings are in line with earlier 
result s 14 ! 42 which report a worse agreement between the 
KS DOS results and the experimental values for stronger 
bound levels. In addition, our results clearly show that 
the agreement between the theoretical and the exper- 
imental spectrum is considerably improved for small 
sodium clusters if the PES is extracted from the true 
excitation energies of the 'daughter' system and not the 
KS DOS. This shows the importance of taking effects be- 
yond the independent-particle picture into account in the 



interpretation of photoelectron spectra. 
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